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ABSTRACT 


Linear  state  signal  feedback  is  used  to  obtain  exponential 
response  from  fourth  order  systems*  Characteristic  equation  roots 
are  selected  to  provide  the  desired  exponential  response  with  con¬ 
straint  on  initial  conditions  and  system  acceleration *  A  digital 
computer  root  locus  program  is  developed  to  determine  feedback 
coefficients  in  a  manner  which  minimizes  the  possibility  of  oscilla¬ 
tory  response  in  the  presence  of  state  sensor  errors*  The  effect 
of  noise  on  the  state  signal  is  investigated  and  a  sample  data  fil¬ 
tering  technique  developed.  A  quasi-optimum  time  technique  utilizing 
second  order  switching  logic  for  initial  control  effort  and  linear 
feedback  for  terminal  control  is  developed* 


ii 


ACKNOWLEDGMENT 


The  author  wishes  to  express  his  appreciation  for  the  assistance 
and  encouragement  given  him  by  Dr*  H.  A,  Titus^  Jroi>  of  the  United 
States  Naval  Postgraduate  School, 


iii 


TABLE  OF  CONTENTS 


Section 

Title 

Pag 

1.0 

Introduction 

i 

2.0 

Linear  Feedback  for  Exponential  Settling 

2 

2.1 

Use  of  the  Uncoupled  Form 

4 

2.2 

Use  of  a  Root  Locus  Study 

8 

2.3 

Consideration  of  Acceleration 

Constraint 

9 

2.4 

Statistical  Considerations 

12 

3.0 

Use  of  Second  Order  Switching 
Quasi-Optimal  Response 

Logic  for 

18 

4.0 

Conclusions 

29 

Bibliography 

30 

APPENDICES 

I. 

Root  Locus  Study 

32 

II. 

Acceleration  Constraint 

39 

III. 

Gaussian  Noise  Generator 

42 

IV. 

Statistical  Study 

45 

V. 

Quasi-Optimal  Control 

51 

iv 


LIST  OF  ILLUSTRATIONS 


Page 


Figure 


1. 

Block  Diagram  of  Aircraft  Longitudinal  Motion 
with  State  Variable  Feedback 

3 

2 . 

Block  Diagram  of  Uncoupled  Form  of  Aircraft 
Longitudinal  Motion  with  State  Variable 

Feedback 

7 

3. 

Block  Diagram  of  General  Fourth  Order  System 
with  Initial  Conditions  on  the  State  Variables 

10 

4. 

Acceleration  of  Fourth  Order  System  with 

Dominant  Characteristic  Equation  Root  of  0,2 

14 

5. 

Time  Response  of  Fourth  Order  System  with 

Noisy  Sensors 

15 

6. 

Time  Response  of  Fourth  Order  System  with 

Noisy  Sensors 

17 

7 . 

Time  Response  for  Saturating  Linear  Control 

23 

8. 

Time  Response  when  Second  Order  Switching  Logic 
is  Applied  to  the  Averaged  Uncoupled  States 

24 

9. 

Time  Response  when  Second  Order  Switching  Logic 
is  Applied  to  the  Uncoupled  State  Fair  Having 
the  Higher  Velocity 

25 

10. 

Time  Response  when  Second  Order  Switching  Logic 
is  Applied  to  the  Uncoupled  State  Fair  Having 
the  Higher  Energy 

26 

11. 

Time  Response  Using  Energy  Priority  and  Revised 
Second  Order  Switching  Logic 

27 

12. 

Time  Response  Using  Energy  Priority  and  Final 

Second  Order  Switching  Logic 

28 

v 


1.0  Introduction 


The  purpose  of  this  study  was  to  investigate  means  of  controlling 
fourth  order  systems  to  obtain  various  types  of  responses.  The  initial 
portion  of  the  study  deals  with  the  use  of  linear  state  signal  feedback 
to  obtain  exponential  settling.  Use  of  the  uncoupled  form  is  investi¬ 
gated  for  use  in  the  feedback  solution  analysis.  The  effect  of  noise 
on  the  state  signals  for  the  linear  feedback  solution  was  investigated 
and  a  sample  data  filtering  technique  developed.  A  root  locus  program 
was  developed  to  determine  feedback  coefficients  in  a  manner  which  will 
minimize  the  possibility  of  oscillatory  response  in  the  presence  of 
state  sensor  errors. 

The  remainder  of  the  study  deals  with  a  quasi-optimum  time 
solution  which  incorporates  switching  logic.  Full  control  effort  is 
used  initially  and  switching  takes  place  at  predetermined  levels  of 
state  variable  combinations.  After  switching,  control  is  allowed  to 
decay  for  a  linear  termination  of  the  solution.  The  objective  here 
is  to  obtain  a  near-  optimum  time  response  with  no  possibility  of  con¬ 
trol  chatter. 

Section  2.0  of  this  report  contains  the  linear  feedback  portion 
of  the  investigation.  Section  3«0  contains  the  quasi-optimum  investi¬ 
gation.  All  synthesis  for  the  study  was  conducted  on  the  school0 s 
Control  Data  Corporation  160k  digital  computer  using  fortran.  All 
programming  will  be  referenced  in  the  text  and  shown  in  the  appendices. 
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2.0  Linear  Feedback  for  Exponential  Settling 


The  response  of  a  system  to  a  set  of  initial  conditions  may  be  con* 
trolled  by  feedback  of  the  system  state  variables.  That  is,  the  forcing 
function  of  the  system  is  made  up  of  predetermined  amounts  of  each  of 
the  system1 s  state  variables.  Thus,  the  characteristic  equation  of  the 
controlled  system  may  be  adjusted  to  obtain  the  desired  response. 

If  an  exponential  response  of  the  system's  position  state  variable 
is  desired,  the  characteristic  equation  will  have  only  negative  real 
roots,  with  a  dominant  root  that  causes  the  desired  exponential  path 
after  the  decay  of  the  non-dominant  roots.  The  proper  feedback  co¬ 
efficients  of  each  of  the  state  variables  may  be  calculated  to  give  the 
desired  characteristic  equation. 

Example 

An  aircraft  landing  flare  is  representative  of  a  class  of  auto¬ 
matic  control  problems  in  which  a  system  has  initial  conditions  of  each 
state  variable  and  it  is  desired  that  the  position  state  variable 
settle  to  zero  in  an  exponential  manner.  In  order  to  design  a  suit¬ 
able  control  system,  a  mathematical  description  of  the  aircraft  longi¬ 
tudinal  motion  is  required.  Assuming  constant  air  speed  and  a  shallow 
glide  angle  leads  to  the  short  period  equation  of  longitudinal  motion 


These  are  written  in  terms  of  the  following  transfer  function 


relating  elevator  position,  d  (radians),  to  aircraft  pitch  angle, 
p  (radians ); 


d(s) 


2 


where 


K  =  short  period  gain 


w  =  short  period  resonant  frequency 
a  =  short  period  damping  factor 
T  =  path  time  constant 

In  order  to  complete  the  mathematical  description  of  the  aircraft, 
the  transfer  function  relating  altitude,  h  (feet),  and  pitch  angle, 
p  (radians),  in  terms  of  velocity  V  (feet  per  second)  and  path  time 
constant  T  is 


h(s)  = 


P(s) 


s(Ts  +  1) 

Combining  (1)  and  (2)  results  in  a  transfer  function  relating 
altitude  and  elevator  position: 


(2) 


h(s)  = 


KV 


5^(s^  +  2as  +  w^) 


d(s) 


(3) 


Letting  the  system  forcing  function  u  be  equal  to  KVd^  a  signal 
flow  diagram  of  the  system  is  shown  in  Figure  1. 


Figure  1.  Block  Diagram  of  Aircraft  Longitudinal  Motion  with 
State  Variable  Feedback 
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The  characteristic  equation  for  the  system  of  Figure  1  is 

sk  +  (b^  +  2a)s^  +  (b^  +  w2)s2  +  b>s  +  b^  ~  0  {h) 

It  is  now  possible  to  solve  for  the  feedback  coefficients 

bf,  b^,  by  and  b^  such  that  the  characteristic  equation  will  have 
the  desired  roots .  Suppose,  for  a  particular  aircraft 
a  =  0*5  and  w  =3  1.0,  and  the  desired  closed  loop  characteristic  equa¬ 
tion  roots  are  s  =-0.l8,  -1„0,  -1.0,  and  -5.0s  giving  the  charac¬ 
teristic  equation 

sk  +  7.18s3  +  12.26s2  +  6.98s  +  0.9  »  0.  (5) 

Equating  (U)  and  (5)  will  then  give  the  desired  feedback  co¬ 
efficients 

bL  =  0.9,  b2  =  6.98,  b3  .  11.26,  b^  »  6.18.  (6) 

This  will  give  a  time  response  solution  for  the  flare,  after  allowing 

a  short  decay  time  (t^  ±  6  seconds)  for  the  non-dominant  roots,  of 
approximately 

h(t)  =  h(tx)  e'°*l8t  (7) 


2 . 1  Use  of  the  Uncoupled  Form 

In  order  to  more  closely  examine  sources  of  instability  in  the 
above  type  of  problem,  and  to  determine  the  proper  variable  for  a  root 
locus  study,  it  is  useful  to  solve  for  the  system*  s  uncoupled  form  |~3 
Rewriting  (3)  in  the  time  domain  gives 


h(t)  +  2  ah(t)  +  w2  h(t)  =  KVd(t) 


(8) 


li 


In  order  to  rewrite  (8)  in  matrix  form  let 

x1(t)  =  h(t) 

=  h(t) 

•  M 

=  h(t) 
x^(t)  =  h”(t ) 
u(t)  =  KVd(t) 

This  leads  to  the  following  matrix  form  for  the  system; 


(9) 


Xg(t) 

*3(t) 


0 

l 

0 

0 

0 

• 

0 

0 

1 

0 

0 

X  = 

0 

0 

°p 

1 

X  + 

0 

0 

0 

2 

-w 

-2a 

1 

(10) 


which  is  defined  as 

x  =  Fx  +  Du  (11) 

The  system  can  now  be  transformed  into  the  uncoupled  form  by 


defining  a  new  variable  y  as 
yL  =  u(s)/s2 
y2  =  u(s)/s 

2  2 
=  u(s)/(s  +  2as  +  w  ) 

2  2 

y^  =  u(s)s/(s  +  2as  +  w  ) 

After  expanding  x  by  partial  fraction  expansion  it  is  noted  that 


(12) 


x  = 


l/w^ 

0 

0 

0 


•2a/w2 

1/w 

0 

0 


(4a2-w2)/w2 

-2a/w 

1 

0 


2a/wp 

-  l/w 
0 
1 


(13) 


which  is  defined  as 
x  = 


(14) 
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Also,  by  solving  for  the  inverse  of  G,  the  expression 


-1 

^  =  G  x 

can  be  written  as 


1  = 


w 

0 

0 

0 


2a  l  o 

w  2a  1 

0  10 

0  0  1 


(15) 


(16) 


In  order  to  draw  a  signal  flow  diagram  of  the  uncoupled  system,, 
(n.)  is  substituted  into  (ll)  giving 
Gj£  =  FG^  +  Du 


(17) 


and  multiplying  by 
1  = 

Solving  (18)  gives 


G-1FG^  +  G-1Du 
the  uncoupled  form 


0 

1 

0 

0 

0 

* 

0 

0 

0 

0 

1 

1  = 

0 

0 

0P 

1 

1  + 

0 

0 

0 

-W2 

-2a 

1 

(18) 


(19) 


which  is  drawn  as  shown  in  Figure  2.  Lines  are  added  to  denote  the 
subsequent  addition  of  state  variable  feedback  loops. 
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Figure  2.  Block  Diagram  of  Uncoupled  Form  of  Aircraft 

Longitudinal  Motion  with  State  Variable  Feedback 

The  characteristic  equation  of  the  uncoupled  system  of  Figure  2 

U  q  2  2 

s  +  (2a+n2+ni+)s3+  (w  +11^+2302+1^)3 


is 

(20) 


2  2 
+  (2an^+w  )s  +  n^w 


=  0. 


Assuming  the  same  aircraft  as  before  with  a  =  0,5  and  w  =  1.0  gives 

s  +  ( l-t-n^+n^  )s  +  ( l+n^+n^+ri^ )  s  +  (n^+n^)s  4-  n^  0  (21) 

Suppose,  prior  to  beginning  the  flare,  the  aircraft  is  descend¬ 
ing  a  -4^  glideslope  with  an  airspeed  of  170  knots  and  that  the 
flare  will  begin  at  100  feet.  Nominal  initial  vertical  velocity  is 
approximately  -20  feet  per  second.  In  order  for  the  glide  angle  to 
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be  approximately  tangent  to  the  flare  path,  a  dominant  characteristic 
equation  root  of  -0.2  is  selected.  That  is,  the  slowest  phase  plane 
eigenvector  is  placed  through  the  nominal  initial  condition  in  the 
x^x^  phase  plane.  Now  the  feedback  coefficients  for  the  control  system 
can  be  calculated.  Suppose  characteristic  equation  roots  are  selected 
at  s  =  -0.2,  -1.0,  -1.0,  and  -5«0.  This  gives  the  following  desired 
characteristic  equation: 

sh  +  7.2s3  +  12.4s2  +  7.2s  +1=0  (22) 

Equating  (21)  and  (22)  gives  the  following  solution  for  the  un¬ 
coupled  state  variable  feedback  coefficients 


n^  =  1.0,  n^  =  6.2,  n^  =  4.2,  =  0 

which  can  be  written  in  matrix  form  as 

N  =  1.0  6.2  4.2  0 

The  solution  for  the  original  state  variable  feedback  coeffi¬ 
cients  can  be  calculated  by  noting 

B  =  NG'1 
which  yields 


{23) 


(24) 


B  = 


1.0  7.2  11.4  6.2 


(25) 


(26) 


2 ,2  Use  of  a  Root  Locus  Study 

In  order  to  become  familiar  with  the  behavior  of  the  roots  of 
the  characteristic  equation  (21),  a  fortran  root  locus  program  was 
written.  Since  the  uncoupled  characteristic  equation  contained 
as  an  element  of  each  of  the  internal  coefficients,  n^  was  used  as 
the  variable  for  the  root  locus.  As  shown  in  Appendix  I,  it  was 
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found  that  (24)  had  all  negative  real  roots  for  the  dominant  root  vary¬ 
ing  from  -0.2  to  -O.366.  Thus  the  previous  solution  is  on  the  boundary 
where  a  slight  error  in  a  feedback  setting  could  cause  two  of  the 
roots  to  leave  the  real  axis  and  result  in  a  oscillatory  component 
in  the  system  response .  For  this  reason^  it  is  better  to  solve  for 
a  dominant  root  of  about  -0.18  and  then  adjust  to  move  the  dom¬ 
inant  root  back  to  -0.2.  Doing  this,  as  shown  in  Appendix  I  gives 


N  = 

[0.9 

5.71 

4.28  0.1  ] 

(27) 

and 

B  = 

[  0.9 

6.61 

10.89  5.81] 

(28) 

These  feedback  coefficients  give  all  real  roots  for  a  dominant 
root  varying  from  -0.18  to  -0.284  when  n^  is  varied  from  4.96  to 
6. 09«  Figure  7  shows  the  system  time  response  resulting  from  use  of 
the  state  variable  feedback  coefficients  of  (28). 

2.3  Consideration  of  Acceleration  Constraint 

The  proper  method  for  selecting  the  non-dominant  root  locations 
of  the  closed  loop  system  characteristic  equation  has  been  ignored 
in  the  previous  discussion.  It  will  now  be  shown  that  system  accelera¬ 
tion  during  exponential  settling  is  a  function  of  both  initial  con¬ 
ditions  and  the  roots  of  the  closed  loop  system  characteristic  eq¬ 
uation.  Thus,  if  the  most  severe  initial  conditions  that  the  system 
can  be  expected  to  be  subjected  to  are  predictable,  an  acceptable 
location  for  the  non-dominant  roots  can  be  determined. 

As  before,  the  dominant  root  should  be  selected  to  place  the 
slow  eigenvector  through  the  nominal  initial  condition  in  the  x^x^ 
phase  plane.  As  the  non-dominant  real  roots  of  the  closed  loop 
system  characteristic  equation  are  moved  to  the  left  from  the  origin 
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in  the  s  plane,  the  acceleration  maximum  during  the  early  portion  of 
the  response  is  increased*  Therefore,  the  non-dominant  roots  can  be 
selected  by  considering  the  maximum  amount  of  acceleration  to  which 
the  system  should  be  subjected*  In  the  aircraft  flare  problem,  the 
proper  acceleration  constraint  would  be  determined  by  structural  con¬ 
siderations  as  well  as  passenger  comfort* 

In  order  to  illustrate  the  effect  of  non-dominant  root  placement, 
a  graphical  solution  of  the  acceleration  on  a  fourth  order  system 
during  exponential  response  was  made*  The  closed  loop  fourth  order 
system  was  generalized  by  lumping  the  feedback  coefficients  with  the 
plant  coefficients  for  each  state  variable*  This  resulted  in  a  gen¬ 
eralized  fourth  order  system  as  shown  in  Figure  3  with  initial  con¬ 
ditions  on  the  state  variables. 


Figure  3*  Block  Diagram  of  General  Fourth  Order  System  with 
Initial  Conditions  on  the  State  Variables 
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Using  signal  flow  techniques,  the  state  variable  transforms  can 
be  expressed  in  terms  of  the  initial  conditions  on  the  state  variables 
As  is  often  the  case  for  mechanical  systems,  the  magnitude  of  maximum 
control  effort  will  be  determined  by  the  acceleration  state  variable. 
The  operational  expression  for  the  acceleration  state  is 


s2h,,(0)  +  h5(0)U3+C|is2)  -  ho(0)(cos+Cl)  -  h1(0)sc1 


(29) 


If  a  purely  exponential  response  on  the  position  state  variable  is 
desired,  the  roots  of  the  closed  loop  system  characteristic  equation 
must  be  real.  This  can  be  expressed  mathematically  as 


(30) 


which  leads  to  the  following  time  expression  for  the  system  acceleration: 


(3D 
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where 


e'Plt 

ei  =  (Pg-Pj  )(p3'-p7^(p'^:p’1T'  etC- 

Curves  may  then  be  plotted  for  various  initial  conditions  and 
characteristic  equation  roots  *  Figure  4  is  a  graphical  solution 
obtained  from  the  fortran  program  in  Appendix  XX.  It  is  a  three 
dimensional  plot  of  system  acceleration  resulting  from  various 
initial  conditions  on  the  position  and  velocity  state  as  well  as 
various  non-dominant  root  locations.  The  dominant  root  location 
is  at  -0,2. 

2  .4  Statistical  Considerations 

This  portion  of  the  study  was  made  to  investigate  the  effect 
of  noisy  state  variable  signals  when  linear  feedback  is  used.  This 
would  be  the  actual  case  in  the  previous  aircraft  flare  example 
since  a  predictable  amount  of  error  due  to  the  aircraft  state  sensors 
could  be  expected.  Radio  altimeter  errors  would  be  caused  by  thermal 
gaussian  noise,  terrain  uneveness,  and  the  effects  of  close  proximity 
to  the  ground  during  the  flare.  The  use  of  an  altimeter  as  the 
prime  sensor  would  necessitate  differentiating  three  times  in  order 
to  obtain  four  state  variables.  Because  of  the  amount  of  noise  that 
could  be  expected  on  the  altimeter  output,  this  would  be  a  highly 
unlikely  approach.  Since  an  accelerometer  would  not  be  sensitive 
to  terrain  uneveness  and  proximity,  a  simulated  accelerometer  output 
was  used  as  the  prime  system  sensor  for  the  computer  synthesis. 

A  simulated  altimeter  output  was  used  as  a  secondary  system  sensor. 
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Figure  4  Acceleration  of  Fourth  Order  System  with  Dominant 
Characteristic  Equation  Root  of  0.2 
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Sensor  noise  was  simulated  using  a  gaussian  noise  generator  [4je 
See  Appendix  III  for  the  mathematical  development  and  computer  check¬ 


out  of  this  noise  generator  * 

A  digital  computer  was  used  to  simulate  the  control  problem*  A 
fourth  order  Runge-Kutta  numerical  integration  algorithm  was  used  to 
simulate  system  motion* 

As  in  the  previous  example,  it  was  desired  that  the  system  be  con¬ 
trolled  over  an  exponential  flare*  The  plant  of  Figure  1  with 
a  as  0*5  and  w  =  1  *0  was  used*  The  feedback  coefficient  of  equation 
(28)  were  selected  to  give  a  dominant  closed  loop  characteristic 
equation  root  at  S  =  -0,2* 

In  the  computer  synthesis,  the  simulated  accelerometer  output^ 
with  gaussian  noise,  is  used  to  calculate  the  four  system  states . 
Smoothing  is  used  to  minimize  the  effects  of  the  added  noise* 
Periodically,  the  altitude  state  is  updated  with  a  simulated  noisy 
altitude  output  which  has  also  been  smoothed*  Updating  is  accomplish¬ 
ed  by  averaging  the  two  altitudes* 

Smoothing  of  both  the  simulated  acceleration  and  altitude  is 
accomplished  by  calculating  a  least  square  error  line  over  a  vari¬ 
able  number  of  past  sensor  outputs*  See  Appendix  IV  for  a  mathe¬ 
matical  development  of  the  smoothing  method* 

Figure  5  shows  the  time  response  of  the  system  in  a  three 
dimensional  graph*  This  multi-curve  graph  shows  the  effects  of 
accelerometer  noise  variation  on  the  response*  Eleven  curves  are 
shown  with  accelerometer  noise  variance  varying  from  zero  to  0*01 
feet  per  second  per  second.  Lines  joining  points  of  equal  time 
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Figure  5  Time  Response  of  Fourth  Order  System  with 
Noisy  Sensors 
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(1  second,  2  seconds,  etc,)  are  shown  to  give  an  illusion  of  depth  to 
the  graph.  For  each  curve  the  altitude  noise  variance  was  0,25  feet. 
Linear  smoothing  over  three  sample  points*  having  a  time  duration  of 
0,025  second  apart,  was  used  for  both  acceleration  and  altitude. 

These  curves  show  that  the  terminal  portion  of  the  response  becomes 
erratic  when  accelerometer  variance  is  0,003  feet  per  second  per 
second  or  more. 

Figure  6  is  similar  to  Figure  5  except  that  the  third  axis  shows 
the  effect  of  using  a  variable  number  of  past  sample  points  for  linear 
smoothing.  Eleven  curves  are  shown  with  the  number  of  sampling 
points  varying  from  zero  (no  smoothing)to  fwelve  ,  Acceleration 
variance  was  0,005  feet  per  second  per  second  and  altitude  variance 
was  0,25  feet.  It  is  seen  that  as  sampling  points  are  increased  to 
five  or  more,  the  data  becomes  sufficiently  stale  so  that  oscillation 
occurs.  As  the  number  of  past  sample  points  is  increased  even  more* 
instability  occurs. 
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Figure  6 


3  *0  Use  of  Second  Order  Switching  Logic  for  Quasi- 


Optimal  Response 

Minimum  time  solutions  for  specific  fourth  order  systems  can  be 
found  in  the  literature  ^5]  -  These  solutions  are  derived  through 
the  use  of  Pontryagin* s  maximum  principle* 

At  this  time,  a  general  solution  for  the  fourth  order  minimum 
time  problem  has  not  been  developed  because  of  the  high  degree  of 
complexity  involved.  This  portion  of  the  investigation  was  conducted 
to  investigate  the  application  of  the  well  developed  second  order 
minimum  time  solutions  to  the  fourth  order  problem  |6j  0 

The  time  optimal  solution  for  fourth  order  systems  requires 
continuous  maximum  control  effort,  For  all  but  discrete  sets  of 
initial  conditions,  three  switching  points  where  control  effort  re¬ 
verses,  are  required *  In  this  study,  a  quasi-optimal  solution  was 
obtained  by  using  switching  logic  for  only  one  switching  point*  The 
philosophy  is  to  use  maximum  control  effort  during  the  initial 
portion  of  the  solution  and  linear  control  for  exponential  settling 
during  the  terminal  portion.  Switching  logic  is  used  to  terminate 
the  maximum  control  effort  at  the  proper  time.  This  switch  must 
take  place  sufficiently  early  to  prevent  excessive  overshoot,  yet 
late  enough  to  result  in  a  reasonably  fast  response.  Control  effort 
is  allowed  to  decay  during  the  terminal  portion  of  the  solution*  This 
linear  portion  then  replaces  the  last  two  switches  of  the  minimum  time 
solution  and  therefor  prevents  chatter  mode* 
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In  section  2.1  it  was  shown  that  the  fourth  order  system  could  be 
represented  in  the  uncoupled  form  as  two  second  order  systems 0  By 
doing  this,  the  well  developed  second  order  switching  lines  of  Titus 
and  Demetry  can  be  utilized  to  give  an  approximate  type  of  switching 
logic  for  the  fourth  order  system. 

Since  two  second  order  systems  result  from  the  uncoupling  of  the 
fourth  order  system,  some  sort  of  priority  must  be  used  to  determine 
which  uncoupled  state  pair  should  be  applied  to  the  second  order  switch¬ 
ing  criteria.  Various  priority  schemes  were  used  including  selection 
of  the  second  order  uncoupled  state  pair  having  higher  velocity,  higher 
energy,  and  variations  of  the  latter.  Also  a  switching  line  which 
averaged  the  states  of  both  the  uncoupled  state  pairs  was  tried. 

The  mathematical  expression  for  minimum  time  second  order  switch¬ 
ing  lines  varies  according  to  the  type  of  plant  eigenvalues  involved 
[«]•  For  this  investigation,  the  second  order  switching  line  for 
null  roots  was  used.  This  switching  line  is 


control 


N  sgn 


xi + 


(33) 


where  N  is  the  saturated  control  effort. 

A  digital  computer  was  used  to  investigate  this  type  of  control. 
See  Appendix  V  for  the  digital  computer  program.  In  each  of  the 
following  cases,  the  plant  of  Figure  1  with  a  -  0.5  and  w  =  1.0  was 
used.  The  terminal  linear  control  was  determined  by  the  feedback 
coefficients  described  in  equation  (28)  which  give  a  dominant  charac¬ 
teristic  equation  root  of  -0.2. 


19 


Figure  7  shows  the  system  time  response  when  only  saturating 
linear  control  is  used.  Initial  system  position  for  each  curve  is 


100  feet.  Initial  system  velocity  is  varied  from  30  feet  per  second 
to  -30  feet  per  second.  Control  saturation  occurs  at  plus  or  minus 
100.  Time  required  for  completion  of  90  percent  of  the  desired  travel 
varies  from  approximately  7  seconds  to  17  seconds  for  the  extreme 
initial  conditions.  These  travel  times  will  be  used  to  evaluate  the 
following  control  systems  which  use  full  control  during  the  initial 
portions  of  each  solution. 

Figure  8  shows  the  system  time  response  when  the  second  order 
switching  logic  of  equation  (33)  is  applied  to  the  average  of  the 
uncoupled  positions  and  the  average  of  the  uncoupled  velocities. 

The  initial  control  equation  therefore  is  changed  to 


control  ~  sgn 


Yl+Y3  + 


Y2  lY2 


+  Y.  Y. 


2N 


(3*0 


Initial  conditions  and  control  saturation  remain  the  same 0  Time  re¬ 
quired  for  completion  of  90  percent  of  the  desired  travel  varies  from 
approximately  3  seconds  to  10  seconds.  This  switching  logic  causes 
the  system  to  respond  approximately  twice  as  fast  as  it  did  when  only 
linear  control  was  used.  However,  this  control  logic  allows  overshoot 
when  the  system  has  negative  initial  velocity.  Therefore,  it  may  be 
unacceptable  for  some  applications. 

Figure  9  shows  the  system  time  response  when  the  second  order 
switching  logic  of  equation  (33)  is  applied  to  the  uncoupled  state  pair 
having  the  higher  velocity.  Initial  conditions  and  control  saturation 
level  remain  the  same.  This  shows  that  the  use  of  the  uncoupled 
velocity  to  determine  the  uncoupled  state  pair  priority  is  not  accept- 
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able.  For  initial  system  velocities  of  10  feet  per  second  or  more,  the 
response  is  identical  to  the  linear  response  of  Figure  7°  For  initial 
system  velocities  of  zero  feet  per  second  or  less*,  excessive  overshoot 
occurs . 

Figure  10  shows  the  system  time  response  when  the  second  order 
switching  logic  of  equation  (33)  is  applied  to  the  uncoupled  state 
pair  having  the  higher  energy.  To  do  this,  each  of  the  uncoupled 
velocity  states  was  normalized  by  dividing  by  system  natural  frequency, 
w.  Then  the  uncoupled  state  pair,  whose  states  were  farther  from  the 
phase  plane  origin,  were  assumed  to  have  the  higher  energy.  Initial 
conditions  and  control  saturation  level  remain  unchanged.  A  fairly 
consistent  overshoot  resulted  from  all  the  initial  conditions  con- 
sidered.  Therefore,  this  control  logic  apparently  has  some  potential. 
However,  the  logic  must  be  altered  to  cause  the  switching  point  to 
occur  sooner. 

Figure  11  is  similar  to  that  of  Figure  10  except  that  the  control 
logic  was  changed  to  cause  the  switching  point  to  occur  earlier.  The 
initial  control  equation  was  changed  to 


control 


“N  sgn 


(35) 


A  considerable  improvement  was  realized  by  changing  the  control  logic. 
The  overshoot  resulting  from  all  initial  conditions  is  now  much 
smaller.  Time  required  for  completion  of  90  percent  of  the  desired 
travel  varies  from  approximately  3  seconds  to  5  seconds  for  the  extreme 
initial  conditions . 
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Figure  12  is  similar  to  Figure  11  except  the  initial  control  logic 
was  again  changed  to  cause  the  switching  point  to  occur  earlier.  The 
initial  control  equation  was  changed  to 


control 


-N  sgn 


y2  Y 


0.8  N 


(36) 


Initial  system  positions  are  100  feet,  70  feet,  and  40  feet.  Initial 
system  velocities  are  30,  0,  and  -30  feet  per  second.  Moderate  over¬ 
shoot  occurs  for  small  initial  position  combined  with  large  negative 
velocity.  Otherwise,  the  system  is  not  extremely  sensitive  to  small 
initial  positions.  This  control  system  would  therefore  be  acceptable 
for  some  applications  where  a  small  amount  of  overshoot  could  be 
accepted  under  extreme  initial  conditions. 
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Figure  7  Time  Response  for  Saturating  Linear 
Control 
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■Figure  8  Time  Response  when  Second  Order  Switching 
Logic  is  Applied  to  the  Averaged  Uncoupled 
States 
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Figure  9  Time  Response  when  Second  Order  Switching 
Logic  is  Applied  to  the  Uncoupled  State 
Pair  Having  the  Higher  Velocity  ' 
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Figure  10  Time  Response  when  Second  Order  Switching 
5  •  Logic  is  Applied  to  the  Uncoupled  State 

pair  Having  the  Higher  Energy 
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Figure  11  Time  Response  Using  Energy  Priority  and 
Revised  Second  Order  Switching  Logic 
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Figure  12  Time  Response  Using  Energy  Priority 

and  Final  Second  Order  Switching  Logic 
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4.0  Conclusions 


By  selecting  the  proper  closed  loop  system  characteristic  equa¬ 
tion  roots,  a  fourth  order  system  may  be  controlled  to  give  a  desired 
exponential  response  within  initial  condition  limitations.  The 
dominant  characteristic  equation  root  will  completely  define  the 
system  output  trajectory  after  decay  of  the  components  associated 
with  the  non-dominant  roots.  Maximum  system  acceleration  is  a 
function  of  initial  conditions  and  characteristic  equation  roots „ 

A  root  locus  study  may  be  used  to  select  characteristic  equation 
roots  which  will  minimize  the  possibility  of  oscillatory  response 
components  in  the  presence  of  state  sensor  errors. 

Within  initial  condition  limitations,  quasi-optimum  time  con¬ 
trol  of  a  fourth  order  system  is  possible  using  second  order  switch¬ 
ing  logic  combined  with  terminal  linear  control. 

It  is  suggested  that  further  studies  be  conducted  to  investigate 
the  use  of  other  types  of  switching  logic  for  optimum  control  of 
fourth  order  system. 
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APPENDIX  I 


ROOT  LOCUS  PROGRAM 


Program  RTLOCUS 
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• . JOB*HARR I S  GRAPH  BOX  228 
PROGRAM  RTLOCUS 

lV(5O^C0NV(50)900'4)'TIPl900’4),OEVt  4’4>  *BIN( 16) , A{ 50) ,J (50) , 

5  FORMAT ( ///2X • 1 4HFEEDBACK  LOOPS, 9X, 23 HCHARAC TER  I  ST  I C  EOUA  TION. 

c£,*  tH2OOT2»  nxt5HROOT3,  1  lX,5HR00T4/2Xi 
R?H^nB7y  AnbocA?43  F?b  aA6X  *  22§!J  f  3X»  24S313X»  2HS2 , 3X,  $HS 1  ,3  X , 

4  IMAG)  ’  °HRE  L  IMAG  REAL  IMAG  R6AL  IMAG  REAL 
PRINT  5 

10  FORMAT ( I5»8F5.3) 

hEis  ord6R0cf1polyn6mul’OEL'*’fsi'p1’ FBLP2lF6l-,,3>F0l-pl* 

A(l-4)  ARE  COEFFS  OF  POLYNOMIAL 
FBLPU-4)  ARE  FEEDBACK  LOOP  COEFFS 

DEL  (  IT4*)  ARE  ^  THE  t  CJHANGE  IN  FEEDBACK  LOOP  COEFFS  FOR  EACH  ITERAT 
tddT^  I  ^NAK,n  -J-Tn !  1  Al*£  R£AL  AND  IMAG  PARTS  OF  UNSORTED  ROOTS 

THE  SORTEDNROOTS  ~4  ARt  ™E  REAL  >ARTS  AND  IMAG  PARTS  0F 
DO  20  K= 1 ,800 
A ( 1 >=1.0 

A(2)=l .0+FBLP2+FBLP4 
A  I  3  1  =  1 .0+FRLP1+FBLP2+FBLP3 
A(4)=FBLP1+FBLP2 
A ( 5 ) =FBLP 1 

CALL  RTPLSUB  I M . A , U , V , CONV 1 
IF  ( K- 1 )  12,12,600 
12  TRP ( K, 1 )=U( 1  ) 

TIP(K, 1 )=V( 1 ) 

TRP(K,2)=U(2) 

TIP(K,2)=V(2) 

TRP ( K , 3 1 =U I  3 ) 

TIP(K,3)=V(3) 

TRP ( K , 4 ) =U  t  4 ) 

TIP (K,4 ) =  VI 4) 

GO  TO  14 

600  DO  605  N=  1 , 4 
DO  605  J=  1 , 4 

l(V(NWlP(K^NiT|^iKilfJ>,‘(U(N,-nP(K-1*J,)  +  (V(N)-TIP(K-l.J)) 
605  CONTINUE  ’ 

. DO  60  1=1 ,4 
I  T=  1 

DO  1000  N=1 ,4 
DO  1 000 J= 1 , 4 
BIN( I T 1 =DEV I N  » J 1 
I  T= I T+ 1 
1000  CONTINUE 
LAG  =  1  • 

'  DO  1020  L=2 » 1 6 

,0,0  ^H^;{76IN,UI  '020.1020.1010 

BIN( 1 ) =B IN( L ) 

B I N ( L ) =TEMP 
L  AG  =  L 

1020  CONTINUE 

30  n2iT°  (30»31»32»33»34,35,36,37,38,39,40,41,42,43,44,45),  LAG 
J=1 

DEVI  1  ,  1  )  =  1 000 . 

DEV 11,21=1000.  • 

DEVI  1 ,31=1000. 

DEVI  1,41  =  1000. 

DEV(2,1 1=1000. 

DEV (3,1 1=1000. 

DEVI4, 1 1=1000. 

_  GO  TO  50 

31  N=  1 
J  =  2 

oevii , i)=iooo. 
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DEV ( 1,2) 
DEV ( 1 ,3) 
DEV ( 1 ,4) 
DEV ( 2 , 2 ) 
DEV ( 3 , 2 ) 
DEV ( 4 , 2 ) 
„„  GO  TO  50 

32  N=  1 
J  =  3 

DEV(1,1) 
DEV( 1 ,2) 
DEV { 1,3) 
DEV( 1 ,4) 
DEV ( 2 , 3 ) 
DEV ( 3, 3 ) 
DEV ( 4 , 3 ) 
GO  TO  50 

33  N=1 
J  =  4 

DEV { 1,1) 
DEV  (1,2)  = 
DEV( 1 ,3) 
DEV  (  1,4)  = 
DEV ( 2 , 4 ) = 
DEV (3,4  )  = 
DEV ( 4 , 4 ) = 
GO  TO  50 

34  N  =  2 
J=1 

DEV ( 2  » 1 )  = 
DEV { 2 , 2 ) = 
DEV (2,3)= 
DEV (2,4)= 
-  DEV (1,1)  = 
DEV ( 3 , 1 ) 5 
DEV (4,1)= 
GO  TO  50 

35  N=2 
J  =  2 

DEV (2,1)= 
.  DEV (2,2)= 
OEV (2,3)= 
DEV ( 2 , 4 ) = 
OEV (1,2)= 
DEV (3,2)= 
DEV ( 4 , 2 ) = 
GO  TO  50 

36  N=2 
J  =  3 

DE V ( 2 , 1  )  = 
DEV (2,2)= 
DEV ( 2 , 3 ) = 
DEV ( 2 , 4 ) = 
OEV (1,3)= 
DEV (3,3)= 
DEV (4,3)= 
„  GO  TO  50 

37  N=2 
J  =  4 

DEV (2,1)= 
DE V ( 2  » 2  )  = 
DEV ( 2 , 3 ) = 
DEV ( 2 , 4  )  = 
DEV ( 1 , 4  )  * 
OEV ( 3 , 4 ) = 
DEV ( 4 , 4 ) » 
GO  TO  50 

38  N=3 
J=1 

DEV (3,1)= 


=1000. 

=1000. 

=1000. 

=1000. 

=1000. 

=1000. 


=1000. 
=1000. 
=1000. 
=1000. 
=1000. 
=1000. 
=1000. 


1000. 

1000. 

1000. 

1000. 

1000. 

1000. 

1000. 


1000. 

1000. 

1000. 

1000. 

1000. 

1000. 

1000. 


1000. 

1000. 

1000. 

1000. 

1000. 

1000. 

1000. 


1000. 

1000. 

1000. 

1000. 

1000. 

1000. 

1000. 


1000. 

1000. 

1000. 

1000. 

1000. 

1000. 

1000. 


1000. 


39 


40 


41 


42 


43 


44 


45 


DEV { 3 ( 2 ) 5 
DEV ( 3 , 3 ) = 
DEVI  3,4)  = 
DEV (1,1)  = 
DEV { 2 , 1  )  = 
DE V ( 4 , 1  )  = 
GO  TO  50 
N=3 
J  =  2 

DE V ( 3 » 1  )  = 
DE V ( 3 , 2 ) = 
DEV (3*3)= 
DEV (3,4  )  = 
DEV (  1,2)  = 
DEV  (2,2)  = 
DEV { 4 , 2  )  = 
GO  TO  50 
N=3 
J=3 

DEV { 3 , 1 ) = 
DE V ( 3, 2  )  = 
DEV { 3, 3 ) = 
DEV{ 3,4)= 
DEV(  1 ,3)  = 
DEV ( 2 , 3 ) = 
DEV (4,3)= 
GO  TO  50 
N=3 
J=4 


1000. 

1000. 

1000. 

1000. 

1000, 

1000, 


1000. 

1000, 

1000. 

1000. 

1000. 

1000. 

1000, 


1000, 

1000, 

1000, 

1000, 

1000, 

1000, 

1000, 


DEV ( 3, 1 ) = 1000 . 
DEV ( 3 » 2 )  =  1 000 . 
DEV{3,3)=1000. 
DEV(3,4)=1000. 
DEV ( 1 ,4)=1000. 
DEV (2,4) =1000. 
DEV ( 4 , 4 ) =1 000 . 
GO  TO  50 
N  =  4 
J=1 

DEV ( 4 , 1 ) = 1 000 . 
DEV { 4  » 2 )  =  1 000 . 
DEV(4,3)=1000. 
DEV{4,4)=1000. 
DEV ( 1 , 1 ) = 1 000 . 
DEV(2,1)=1000. 
DEV ( 3 , 1 ) = 1 000 . 
GO  TO  50 
N  =  4 
J  =  2 


DEV (4,1) 
DEV (4,2) 
DEV (4,3) 
DEV (4,4) 
DEV{ 1 ,2) 
DEV (2,2) 
DEV (3,2) 
GO  TO  50 
N  =  4 
J  =  3 

DE V ( 4 , 1 ) 
DEV { 4 , 2 ) 
DEV  (4,3) 
DEV ( 4 , 4 ) 
DEV( 1 ,3) 
DEV (2,3) 
DEV (3,3) 
GO  TO  50 
N=4 
J=4 

DEV (4,1) 


>1000, 

>1000, 

^1000, 

1000, 

1000, 

1000, 

1000, 


:1000, 

1000. 

1000, 

1000, 

1000, 

1000, 

1000, 


=1000, 


DE V ( 4 , 2 ) = 1 000 . 

DEV(4,3)=1000. 

DEV ( 4 , 4 )  =  1 000 . 

DEV ( 1 , 4 ) = 1 000 . 

DEV(2,4)=1000. 

DEV ( 3 , 4 ) = 1 000 • 

GO  TO  50 

50  TRP(K,J)*U(N) 

TIP(K,  J)=V(.M) 

60  CONTINUE 

14  CONTINUE 

15  FORMAT! 1 X , 4F5 . 2 , 3X, 5F5 . 2 , 4X , 8F8. 3 ) 

PRINT  15,FBLP1,FBLP2,FBLP3,FBLP4,A(1 ),A(2),A(3),A(4),A(5), 

1  T  R  P  ( K ,  1  )  » T I P.  (  K ,  1  ) ,TRP(K, 2),TIP(K,2),  TRP  (  K,  3  )  ,  TI  P  (  K ,  3  )  , 

2TRP(K,4) ,TIP(K,4> 

FBLP1=FBLP1 +DEL1 
FBLP2=FBLP2+DEL2 
FRLP3=FBLP3+DEL3 
FBLP4=FBLP4+DEL4 
NUMPTS=K 

20  CONTINUE 

00  22  L= 1 i M 
00  21  K= 1 , NUMPTS 
TRP ( K ) =TRP ( K , L ) 

T  I  P ( K )  =  T I P ( K, L ) 

21  CONTINUE 

CALL  GRAPH  ( NUMPTS , TRP , T I P, 8 ) 

22  CONTINUE 
END 

SUBROUTINE  RT PLSUB ( N, A, U , V, CONV ) 

DIMENSION  A (50) ,H(50> , B( 50) ,C(50),D(  50) , E(50) ,U(50> , V(50 ) ,CONV( 
F= 1 0. 0 
L  =  25 


54 

52 

100 


101 


150 

151 


152 

200 

205 

201 


I F ( N )  52,54,52 
PAUSE  30 
NP3=N+3 
B ( 2  )  =0«  0 
B( 1 )=0.0 
C  (  2  )  =  0. 0 
C(  1  )=0.0 
D(2)=0.0 
E ( 2 ) =0 . 0 
H(2)=0.0 
DO  101  J=3,NP3 
H ( J ) = A ( J-2 ) 

T=1 .0 

SK= 1 0 • 0**F 

I F ( H ( NP3 ) )  200,151,200 
U ( NP3 ) =0. 0 
V ( NP3 ) =0. 0 
CONV ( NP3 ) =SK 
NP3=NP3- 1 

IF(NP3) 152, 152, 150 
STOP  152 

IF(NP3-3)205,51 ,201 

STOP  205 

PS=0.0 

QS*0 . 0 

PT=0 . 0 

QT=0. 0 

S=0.0 

REV=1 .0 

SK= 1 0. 0**F 


I 


IFINP3- 4)206, 202,203  • 

206  STOP  206 

202  R=-H(4)/H(3) 

GO  TO  500 

203  DO  207  J  =  3»  NP3 
IF(H(J) 1204,207,204 

204  S=$+LOGF ( ABSF ( H ( J ) ) ) 
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207 


208 

210 

250 


251 

252 

253 


254 

255 


256 

257 

258 
300 

350 

351 

352 

353 
356 

354 

400 

401 
403 

450 

451 

452 
460 

453 


454 


455 

456 

457 

458 
490 


491 

492 

500 

501 


CONTINUE 
FPN 1 =N+ 1 
S=EXPF( S/FPNl ) 

DO  208 J=3 , NP3 
H(J)=H(J)/S 

IF(A8SF(H(4)/H(3) ) -A8SF ( H( NP3- 1 ) /H(MP3) ) 1250,252,252 
T  =  -  T 

M=(NP3-4)/2  +  3 
DO  251  J=3,M 
S=H { J 1 
JJ=NP3-J+3 
H ( J } =H( J J ) 

H  (  J  J  )  =  S 

IF(QS)  253,254,253 

P  =  PS 

Q=QS 

GO  TO  300 
HH2=H { NP3-2 ) 

I F ( HH2 )  256,255,256 

Q=1 .0 

P=-2.0 

GO  TO  257 

0=H(NP3)/HH2 

P= ( H ( NP3-1 ) -Q*H ( NP3-3 ) 1/HH2 
IF (NP3-5)258,550,258 
R  =  0 . 0 

DO  490  1=1, L 

DO  351  J=3, NP3 

8{ J)=H( J)-P*B(J-1 )-Q*B( J-2) 

C( J)=B( J ) -P*  C { J- 1 )-Q*C( J-2) 

I F { H ( NP3- 1 ) ) 352,400,352 
IF<  BINP3-1 1 1353,400,353 
AVH8l=ABSF(H{NP3-l )/B(NP3-l) 1 
IF(AVHBl-SK)450,354,354 
B(NP3)=H(NP3)-Q*B(NP3-2) 

I F ( B ( NP3 1 1401,550,401 
AVHB2=ABSF (H(NP3)/B(NP3)  1 
IF(SK-AVHB2)5 50 *450,450 
DO  451  J=3, NP3 
D( J)=H(J1+R*D( J-l 1 
E( J)=D( J)+R*E( J-l ) 

I F ( D ( NP3 1 1452,500,452 
AVHD3=ABSF(H(NP3)/D(NP3) ) 

I F( SK-AVHD3 1500,453,453 

CC2=C ( NP3-2 1 

CC3=C ( NP3-3 1 

C ( NP3- 1 )=-P*CC2-Q*CC3 

CC 1 =C ( NP3-1 1 

S=CC2*CC2-CC1*CC3 

IF(S)455,454,455 

P=P-2 . 0 

Q=Q* ( Q+ 1 .01 

GO  TO  456 

P=P+  (  B  (  NP3-1  )*CC2-B(NP3)-»CC3)/S 
Q=Q+ ( -B ( NP3- 1 1*CC1+B(NP3)*CC2)/S 
I F ( E ( NP3- 1 1 1458,457,458 
R=R- 1.0 
GO  TO  490 

R=R-D(NP3)/E(NP3-1 1 
CONTINUE 


PT  =  P 
QT=Q 

I F ( RE V ) 49 1 , 492 , 492 
SK=SK/ 1 0 . 0 
REV=-REV 
GO  TO  250 
I F ( T 1 50 1 ,502,502 
R= 1 . 0/R 
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ojcvjro 


502 


503 

550 

551 

552 

560 

553 


554 


555 

556 

557 


561 

558 

51 


NP=NP3-3 
l){  NP)=R 
V ( NP ) =0. 0 
CONV ( NP ) =SK 
NP3=NP3- 1 
DO  503  J  =  3,  NP3 
H  (  J ) =D { J ) 

IF(NP3-3>300,51 ,300 

IF(T) 55 1*552,552 

P=  P/Q 

Q= 1 . 0/Q 

PP2=P/2. 0 

QMPSQ=Q-PP2*PP2 

IF(0MPSQ)554,554,553 

NP=NP3-3 

U { NP ) =  -PP2 

U { NP- 1 ) =  —  P  P  2 

S=SQRTF ( QMPS  Q ) 

V  (  NP ) =S 
V { NP-1 ) =  -S 
GO  TO  561 
S=SQRTF{-QMPSQ) 
NP=NP3-3 

I F ( P ) 555 ,556*556 
U(NP)=-PP2+S 
GO  TO  557 
U { NP ) =-PP2-S 
U ( NP-1 ) =Q/U ( NP ) 

V ( NP ) =0 . 0 
V( NP-1 )=0.0 
CONV { NP ) =SK 
CONV (NP-1 ) =SK 
NP3=NP3-2 
DO  558  J=3, NP3 
H  (  J  )  =B ( J ) 

GO  TO  200 

RETURN 

END 

END 
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APPENDIX  II 


Graphical  Solution  of  System  Acceleration 


Program  SETROOT 
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..JOB*HARRIS  BOX  228 
PROGRAM  SETROOT 

DIMENSION  ACC  1 (900)  , ACC2 ( 900 ) , ACC3 ( 9 00 ) , ACC 4 ( 900 ) , ACC  TOT (900) , 

1 T ( 900 ) , X ( 900 ) , Y ( 900 ) 

DO  200  K= 1 . 30 
5  FORMAT  (9F10.5) 

READ  5,ALTIC,SINKIC, ACC  I  C ,  BERK  IC  ,  P  1,  P2.P3,P4 

C0EF1 =P1 *P1  *  ( ( ALT IC*P2*P3*P4)+SlNKi:  * (P2*P3  +  P2*P4  +  P3*P4) +ACCIC* 
1 ( P2+P3+P4 )+BERK I C ) / ( ( P2-P1 )* ( P3-P1 )* ( P4-P1  )  ) 

C0EF2=P2*P2*  (  (ALTIC*P1*P3*P4)+SINKI3*(P1*P3+P1*P4+P3*P4)  +ACCIC* 
1  (Pl  +  P3  +  P4)  +  BERKIC)/(  (  P  1  —  P  2  )  ■*  (  P3-P2  )*  (P4-P2)  ) 

C0EF3=P3*P3*( ( ALT  I C*P  1 *P 2*P4 )  +  S I NK i; * ( Pi *P2  +  P 1  * P4  +  P2* P4) +ACCIC* 
1  (Pl+P2+P4)+BERKIC)/(  (P1-P3)*(P2-P3)*  (  P4-P3)  ) 

C0EF4=P4*P4* ( (ALTIC*Pl«P2*Pi)+SINKi:*(Pl*P2+P1*P3+P2*P3) +ACCIC* 
1  (P1+P2  +  P3)+BERKIC)/ ( ( P1-P4 )*( P2-P4)t (P3-P4)  ) 

T ( 1 )=0.0 

DO  100  J  =  1 »  300 

ACC  1 (J)=C0EF1*EXPF(-P1«T(J)  ) 

ACC2( J)=C0EF2*EXPF(-P2*T ( J)  ) 

ACC3(  J)=C0EF3*EX-PF(-P3*T(J)  ) 

ACC4(J)=C0EF4*EXPF(-P4*T(J  )  ) 

ACCT0T( J)=ACC1 ( J)+ACC2( J )+ACC3( J)+A:C4( J ) 

X(J)=T(J)-SINKIC *0.5- 10.0 
Y( J)=ACCTOT( J )-S INK IC*0. 5*5. 0-50.0 
T( J+l )=T( JJ+0.025 
100  CONTINUE 

CALL  GRAPH  (200fX,Y,8) 

200  CONTINUE 
END 
END 
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APPENDIX  III 


Gaussian  Noise  Generator 


Algorithm 

Given  a  uniform  distribution  p(x)  =  -  „5;Sx  <5 

E(x)  as  U  =*  0.0 

Var(x)  =  \  (x=u)^p(x)dx  =  l/l2 

A  5 

«r  12 

Let  R  as  c  2L.  x 
"  1 

e(r)  -  0.0 

?  «C  12 

Var(R)  =  Var(cnx)  =  c  Var(  x) 

P  1 

=  c  n  Var(x) 

2 


and 


p(R)  ^  Normal(0>c  )  by  the  Central  Limit  Theorem 


Computer  Program  Flow  Chart 


enter 


/ sum  12  rand, 
/numbers  having 
(  unif.  dist. 

\  from  “0.5 
\  to  0.5 

multiply 
by  desired 
deviation 


42 


$ 


% 


Graph  1  Distribution  of  1000  Random  Numbers  versus  Theoretical 
Distribution  with  Standard  Deviation  of  Unity 


i+3 


• • JOB«HARR I S  BOX  228 
PROGRAM  CHKRANO 

DIMENSION  RAND(5000),XMAG(500),BIN(300),XLINE(500) ,YLINE (500) 

VA  R  =  4 . 0 
START=71 .3921 
NUMPTS=5000 
NUMB  I NS= 1 00 

XNUMPTS=NUMPTS  ‘  l 

C  CALL  RANDOM  NUMBER  GENERATER  AND  PLACE  IN  ARRAY 

DO  10  1=1, NUMPTS 

CALL  RANDOM  { RAND(I) , START, VAR ) 

RANTOT=RANTOT  +  RA'ND(  I  )  t 

10  CONTINUE  i 

C  CALCULATE  MEAN  AND  DEVIATION 

DEV=0 . 0 
TOT=NUMPTS 
RANME AN= RANT OT/ TOT 
DO  20  1  =  1 » NUMPTS 

20  DEV  =  DEV+( RANMEAN-RANDt I )  )*( RANME AN-R AND (  I  )  ) 

DEV=SQRTF( DEV/TOT) 

30  FORMAT { / / » 6H  ME AN=  F10.6,5X,11H  DEVIATION=  F10.6) 

PRINT  30 , RANME AN  f  DEV 

C  SORTS  THE  ORDERED  RANDOM  NUMBERS  INTO  240  BINS,  INCREMENTS  BIN! 

C  FOR  EACH  NUMBER  THEREIN 

XMAG ( 1 )=-3.0«SQRTF(VAR) 

UMBINS=NUMBINS 
BININC=~2.0«XMAG( 1 ) /UMBINS 
DO  60  L= 1 , NUMB  I  NS 
XMAGIL+l  )=XMAG(  D+BININC 
B I N { L ) =0 . 0 

DO  50  1=1. NUMPTS  '  • 

IF  (RANO( I)-XMAG(L) )  50,45,45 

45  IF  ( RAND (I)-XMAG(L+1))  47,47,50 

47  BIN(L)=BIN( L  )  +1 • 0 / ( BININC*XNUMPTS) 

50  CONTINUE 
60  CONTINUE 

CALL  GRAPH  ( 240 , XM AG , B I N , 8 ) 

C  COMPUTE  ACTUAL  NORMAL  CURVE 

PI=3. 141593 
XDEV=SQRTF( VAR) 

ACT0R=SQRTF(2.0*PI ) 

C0EF=1 .0/ ( ACTOR*XDEV) 

XLINEt 1 ) =-3 . 0  *XDEV 
DO  70  L= 1 ,300 

XLINEIL+l )=XLINE(L )+XDEV/50.0 

YLINE(L)=CDEF*EXPF(-O.5*(XLINE(L)/X0EV)*(XLINE(L)/XDEV)) 

•  70  CONTINUE 

CALL  GRAPH  ( 300 , XL  I NE , YL I NE , 8 ) 

END 

SUBROUTINE  RANDOM  ( QUANT, START, VAR) 

C  GENERATES  RANDOM  NUMBERS  WITH  NORMA.  DISTRIBUTION 

TOT  ALCL  =  0 . 0 
FACT0R=12. 45321 
DO  91  J  =  1 ,  12 

XNUM= (ST ART +6 1.9) *F ACTOR 
NUMX=  XNUM 
X I  NT  =  NUMX 
START=XNUM-XINT 
TOTALCL=TOTALCL+( START-0.5 ) 

91  CONTINUE 

QUANT=TOTALCL*SORTF ( VAR) 

RETURN  . 

END 

END 

0  8  8 
2  NOISE  DISTRIBUTION 
2  HARRIS  BOX  228 

1 

3 


hk 


APPENDIX  IV 


Linear  Smoothing  of  Noisy  Data 


1 .  Mathematical  Development 

...  <rn  ,  \2  ,2 

minimize:  /  (y  )  =  a 

L-  y  werror' 

I 

define  line  y  -  Ax  -  B 

n 


then  d  =  ^  (y,.  -  Ax^.  -  B)^ 

(1)  dd~  /dA  =  2_  2(y.  “  Axi  "  B)(-x..) 

1 

(2)  dd2  /dB  -  ^  2(y«  -  Ax<  -  b)(-U 

1 

y n  5n 

y±  =  ny>  (_  x.  =  nx 


0 


0 


(2)  B  = 


(1)  A  = 


y  -  Ax 
n 


i 


x.y.  -nyx 


;>nx2 

1  1 


-nxx 
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Program  NXFLARE 
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..JOB  HARRIS  BOX  228 
PROGRAM  NXFL ARE 

THIS  PROGRAM  IS  TO  INVEST IGATE  LINEAR  CONTROL  WITH  NOISY  DATA. 
INITIAL  CONDITIONS  ARE  ASSUMED  KNOWN .  ONLY  ALTITUDE  AND  ACCELE 
SENSORS  ARE  USED.  CHARACTERISTIC  EDUAT  1  ON- ROOTS  ARE  .2,2,2 ,2- 

DIMENSION  X(  4 ) ,XP( 4.900) , TP! 900) ,XJ(  900) , YJ ( 900) ,XN( 4, 900) , 

1  XFR0M3! 4,900 ) .XFROMl  (  900  ) ,  X1N0ISE!  9D0  ) ,  X3N0I  SE  (  900) ,  XRU ) , 

2XX ( 20 , 1 1 ) ,YY(20, 11 ) 

START=71 .3921 

5  READ  3,N,NSAMP,X1VAR, X3VAR,T0,TF,DT,  (X( J  ),J=1,N) 

SAMP  =  NS AMP 

N=ORDER,  T0= I N I T I AL  TIME,  TF=F INAL  TIME,  DT=TIME  INC,  X=STATES 
T  =  TO 

3  FORMAT  (2I5,2F10.5/8F10.5) 

IFIN-13)  6,30,30 

6  PRINT  200, T, ( X C J),J  =  1,N) 

L  =  0 
M=M+  1 
NC0UNT=40 
11  =  1 

TP ( I  1 ) =T 
DO  2  J= 1 , N 
XR ( J ) =X ( J ) 

2  XP ( J , II ) =X ( J  ) 

GO  TO  20 
309  11=11+1 
TP ( II ) =T 
DO  12  J=  1  ,  N 
12  XP( J,I1 ) =X( J  ) 

CALL  RANDOM  ( X3N0 ISE (  I  1  )  , START, X3V AT ) 

XN3  =  X(3)+X3N0ISE(  I  1  ) 

CALL  RANDOM  (XlNOISE( 11 ), START, XlVAR) 

XN1  =X ( 1 )  +  Xl  NOISE! II ) 

IF  ( T-S AMP*OT )  50,60,60 

50  XR ( 1  )=X( 1 ) 

XR ( 2 ) =X ( 2 ) 

XR ( 3 ) =X ( 3 ) 

XR ( 4 ) =X ( 4 ) 

XN ( 3 , 1  1 ) =XN3 
XN ( 1  , 1 1- )  =XN  1 

C  NEXT  TWO  LINES  SET  UP  INITIAL  VALUES  FOR  INTEGRATION  IN  SMOOTH  SUBR 

XFR0M3 (2,11 >  =X( 2 ) 

XFR0M3 ( 1 , 1 1 ) =X ( 1 ) 

XNl =0.0 
XN3=0 . 0 
GO  TO  310 
60  XN ( 3 , 1 1 ) =  X N 3 

XN ( 1 , I  1 ) =XN 1  . 

XN 1 =0 . 0 
XN3  =  0 . 0 

CALL  SMOOTH  ( NSAMP , DT , 3, II , TP , XN , XFT 0M3 ) 

JC0UNT  =  J COUNT  + 1 
IF! 10-JC0UNT)  80,80,70 
70  XR! 1 )=XFR0M3( 1,11) 

XR ( 2 ) =XFR0M3 ( 2 , II ) 

XR ( 3  >  =XFR0M3 (3,11) 

XR(4)=XFR0M3(4, II ) 

GO  TO  310 
80  JCOUNT  =  0 

C  USE  SMOOTH  ONLY  EVERY  TENTH  ITERATION  FDR  XI  THEN  TAKE  AVERAGE  POSI 
C  FOR  NEXT  ITERATION 

CALL  SMOOTH  { NS AMP , DT , 1 ,  1 1 , TP , XN , XOJ T ) 

XFROMl (II ) =XOUT 


kT 


11  CONTINUE 

1 J=1  ,4)  ,  XFROM 1 (II) 

XR  (  1  )  =  (  XFROM3  (  1 ,  I  D  +XFROM1  (  11  )  >/2. 

XR ( 2 ) =XFROM3 ( 2  » I  1 ) 

XR ( 3 ) =XFROM3 ( 3 , 1 1 ) 

XR ( 4 ) =XFROM3 ( 4, I  1 ) 

GO  TO  310 

310  IF  ( TF.-T  )  10,20,20 

200  FORMAT  ( /3X,4HTIME,6X,4HX( 1 ) ,6X,4HX(  2 ) , 6X , 4HX ( 3 ) , 6X , 4HX( 4) ,6X, 
17HX1 NOISY, 3X, 7HX3NO ISY ,3X,7HX1FR0M3,  3X, 7HX2FR0M3, 3X, 7HX3 FR0M3 , 
23X, 7HX4FR0M3,3X,7HX1 FROM  1 // 12F1 0.5)  ' 

100  FORMAT  (12F10.5) 

10  00  15  1=1,11 

YJ(  I  )=XP{ 1 , I ) +X3VAR*4000. 

XJ(  I  ) =TP ( I ) +X3VAR« 1 000 . 

IF  ( NCOUNT-40 )  14,160,160 

160  L=L  + 1 

YY ( L , M ) =YJ ( I ) 

XX ( L , M ) =X J ( I ) 

NC0UNT=0 

14  NC0UNT  =  NC0UNT  + 1 

15  CONTINUE 

CALL  GRAPH  (I1,XJ,YJ,8) 

GO  TO  5 

30  DO  35  L= 1  , 2 1 
DO  32  M=1, 1 1 

YJ ( M ) =YY ( L, M ) 

X J ( M ) =XX ( L » M ) 

32  CONTINUE 

CALL  GRAPH  ( 1 1,XJ,YJ,8)  • 

35  CONTINUE 
STOP 

20  CALL  RKUTTA  ( N, T , XR , X , DT ) 

T =T  +  DT 
GO  TO  309 
END 

SUBROUTINE  RANDOM  ( QUANT , START , VAR ) 

C  GENERATES  RANDOM  NUMBERS  WITH  GAUSSIAN  DISTRIBUTION 

TOT  ALCL=0 • 0 
FACT0R=12. 45321 
DO  91  J=1  , 12 

XNUM= (START+61.9) * FACTOR 

NUMX=XNUM 

X  I  NT  =  NUMX 

ST  ART=XNUM-X I  NT 

TOTALCL=TOTALCL+(START-0.5) 

91  CONTINUE 

QUANT=TOTALCL*SQRTF ( VAR ) 

RETURN 

END 

SUBROUTINE  RKUTTA  ( N, T , XR, X , DT ) 

DIMENSION  X(4),AK(4,4),XD0T(4),XC(4) , C ( 4 ) , XR ( 4 ) 

C( 1 )=0.0 
C ( 2 ) =0.5 
C ( 3 ) =0.5 
C  (  4  )  =  1 . 0 
DO  4  1=1,4 
TC=T+C ( I ) *DT 
DO  2  J= 1 , N 

2  XC(J)=XR(J)+C( I)*AK(I-1  ,J) 

CALL  DERI  V  (XC,XDOT) 

DO  4  J= 1 , N 
4  AK ( I, J)=DT*XDOT(  J) 

DO  3  J=1  ,N 

3  X(J)=XR(J)+(AK( 1,J)+2.*AK(2,J)+2.*A<(3,J)+AK(4,J) )/6. 

END 

SUBROUTINE  DERI V  (X,XDOT) 

DIMENSION  X ( 4  )  , XDOT ( 4 ) 

XDOT ( 4 ) =-l . 0*X ( 3 ) - 1 .  0*X(  4J+C0NTR0L 
XDOT ( 3 ) =X { 4 ) 
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XDO  r ( 2 ) =X ( 3 ) 

XDOTM  )  =X  (  2  ) 

CONTROL=-0.9*X( 1 )  -6.61*X(2)  -10.89*X(3)  -5.81*X(4) 

END 

SUBROUTINE  SMOOTH  (  NX , OX  ,  J  ,  1 1  ,  X, Y,  X<  ) 

DIMENSION  X ( 900 ) , Y ( 4 , 900 ) , X X ( 4 , 900 ) , A ( 3, 4 ) , ANSWER ( 4 ) 

DO  2  J=  1 , 3 
DO  2  1=1 ,4 
A ( J , I ) =0 • 

2  CONTINUE 
FN=NX 

DO  5  1=1, NX 
I N  =  I  —  1 
I T- 1 1  —  IN 

A(3,4 )=A(3,4)+Y ( J, IT) 

A!2,4)=A(2,4)+Y( J, IT)*X< IT) 

A( 1 »  4 ) -A ( 1»4)+Y(J,IT)*X(  IT ) *  X ( IT) 

A ( 2 , 3 )  =  A ( 2»  3 ) +X { IT) 

A(2»2)=A(2,2)+X! I T ) *X ( IT  ) 

A( 1 ,2)=A( 1 ,2)+XI  IT)*X< IT  )*X(  IT) 

A( 1  ,1  )  =  A ( 1, 1  )+X{ I T ) *X (  IT ) *X ( I T  ) *X( IT ) 

5  CONTINUE 

All , 3  )  =A { 2 , 2 ) 

A(2, 1  )=A( 1 ,2) 

A(3, 1 )=A(2,2) 

A  (  3 , 2  )  =  A  (  2 , 3  ) 

A ( 3 , 3  )  =FN 
NUMBER=3 

CALL  J0RDAN2  < A, NUMBER, ANSWER ) 

IF  ( J-2 )  8,8,10 

10  XX ( 3, II ) = ANSWER! 1 )*X( I  1 ) *X (II ) +ANS WE  R ( 2 ) *X ( I  1 )+ ANSWER! 3) 
XX! 4, II ) =2.* ANSWER! 1)*X(  II )+ANSWER!2) 

XX(2, II )=XX( 2, 1 1-1 )+XX! 3, 1 1 )*DX 
XX!1,I1)=XX(  1  ,I1-1)+XX(2,I1 
RETURN 

8  XX=ANSWER ( 1 )*X(I1)*X(I1) +ANSWER! 2 ) *< (ID  +ANSWER(3) 

END 

END 
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APPENDIX  V 


Use  of  Second  Order  Switching  Logic 
for  Quasi-Optimum  Control 


Program  OPTFOUR 
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• . JOB  HARRIS  BOX  228 
PROGRAM  OPTFOUR 

DIMENSION  X(4)  ,  X.P  (  4 , 900)  ,  TP  (900)  ,XJ[  900) ,YJ(900) 

5  READ  3,N,T0,TF,DT,(X(J),J=1,N) 

T=T  0 

3  FORMAT  (  I5»/8F10.5) 

IF  (N-13)  6,30,30 

6  CONTINUE 
11  =  1 

TP( II )=T 
DO  2  J=  1  ,  N 
2  XP( J, II )=X( J) 

GO  TO  20 

309  11=11+1 
TP( II ) =T 

DO  12  J=1 ,N 
12  XP( J,I1 ) =X( J ) 

310  IF  (TF-T)  10,20,20 
10  DO  15  1=1,11 

YJ( I )  =  XP(  1  ,1  ) 

XJ( I ) =TP (  I) 

15  CONTINUE 

CALL  GRAPH  (II ,XJ,YJ,0) 

GO  TO  5 
30  STOP 

20  CALL  RKUTTA  (N,T,X,DT) 

T=T  +  DT 
GO  TO  309 
END 

SUBROUTINE  RKUTTA  (N,T,X,DT) 

DIMENSION  X(4),AK(4,4), X DOT (4),XC(4)  ,C(4) 

C( 1 )=0.0 
C ( 2 ) =0. 5 
C ( 3) =0.5 
C  (  4  )  =  1 . 0 
DO  4  1  =  1  ,4 
TC=T+C(I)*DT 
DO  2  J  =  1 »  N 

2  XC( J)=X( J)+C( I ) *AK( 1-1 , J  ) 

CALL  DERIV  (XC,XOOT,T) 

DO  4  J= 1 , N 

4  AK ( I , J ) =DT *XDOT ( J ) 

DO  3  J=1 ,N 

3  X ( J ) =X ( J ) + ( AK( 1 ,J ) +2 . * AK ( 2, J )  +2 .  »AK( 3 , J ) +AK ( 4 , J ) ) /6* 

END 

SUBROUTINE  DERIV  (X,XDOT,T) 

DIMENSION  X ( 4 ) , XDOT ( 4 ) 

C  SELECTS  UNCOUPLED  STATE  PAIRS  HAVING  HIGHER  ENERGY 

C  ITITIALIZE  FOR  NEW  CURVE 

IF  (T)  2,2,4 
2  JFL AG=0 
KFL AG=0 

4  XD0T(4)=-1 .0*X( 3 ) - 1 .0*X( 4 ) +CONTROL 
XDOT { 3 ) =X ( 4 ) 

XDOT ( 2) =X ( 3) 

XDOT ( 1 )=X(2) 

SAT= 100. 

C  IF  SWITCHING  POINT  PREVIOUSLY  REACHED  USE  LINEAR  CONTROu 

IF  ( JFLAG )  6,6,32 

C  SWITCHING  POINT  NOT  PREVIOUSLY  REACHED 

C  CALCULATE  UNCOUPLED  STATE  VARIABLES 

6  Y 1 =  X ( 1 )  +  X ( 2 ) +X ( 3 ) 

Y2=X(2)+X(3)+X(4) 

Y3=X ( 3 ) 

Y4=X ( 4 ) 

SOMEG= 1 . 

EN12=S0MEG*Y1*Y1+Y2*Y2 
EN34=S0MEG*Y3»Y3+Y4»Y4 
IF  (ABSF(EN12)-ABSF(EN34)  )  10,15,15 

10  Y2=Y4  * 
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OCM 


Y 1  =  Y  3 

15  CONTINUE 

C  CALCULATE  OPTIMUM  CONTROL 

C0N1R0L  =  -SICNF( SAT ,Y1 +Y2*ABSF(  Y2 )  / ( 0 . 8*SAT ) ) 

C  SENSE  WHEN  OPTIMUM  CONTROL  CHANGES  SIGN 

IF  (CONTROL)  50,30,30 

C  CONTROL  POSITIVE  SET  J  FLAG  FOR  LINEAR 

30  JFL AG= 1 

100  FORMAT  ( 1  OH  SWITCH  AT , F 1 0 . 5 , 4HFE ET ) 

PRINT  100, X {  1  ) 

C  CALCULATE  LINEAR  CONTROL 

.  32  CONTROL=-0.9*X(  l')-6.61*X(2)-10.89*X{  3)— 5.'81*X(4) 

C  CHECK  FOR  SATURATION 

IF  (CONTROL-SAT)  40,40,35 
35  CONTROL=SIGNF(SAT, CONTROL) 

40  KFLAG=KFLAG+ 1 

IF  (KFLAG-1)  45,45,50 
45  PRINT  200, X( 1 ) 

200  FORMAT  (  1  5H  UNSATURATEO  AT , F 1 0 .5 , 4HF EET ) 

50  RETURN 
ENO 
END 
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SUBROUTINE  DERIV  (X,XDOT,T) 

DIMENSION  X ( 4 ) ,XD0T(4 ) 

C  SELECTS  UNCOUPLED  STATE  PAIRS  HAVINS  HIGHER  .VELOC I TY 

C  ITITIALIZE  FOR  NEW  CURVE 

IF  (T)  2,2,4 
2-  JFL AG  =  0 
KFLAG=0 

4  XDOT(4)=-1.0*X(3)-1.0*X(4)+CONTROL 
XOOT ( 3 ) =  X ( 4 ) 

XOOT ( 2 ) =X ( 3 ) 

XDOT( 1 ) =X ( 2 ) 

SAT  =  100. 

C  IF  SWITCHING  POINT  PREVIOUSLY  REACH; 0  USE  LINEAR  CONTROL 

IF  ( JFLAG )  6,6,32 

C  SWITCHING  POINT  NOT  PREVIOUSLY  REACHED 

C  CALCULATE  UNCOUPLED  STATE  VARIABLES 

6  Yl=X( 1 )+X(2)+X(3) 

Y2=X(2)+X(3)+X(4) 

Y3=X ( 3 ) 

Y4=X (4) 

IF  IABSFI Y2)-ABSF(Y4) )  10,15,15 
10  Y2=Y4 
Y 1  =  Y3 

. 15  CONTINUE 

C  CALCULATE  OPTIMUM  CONTROL 

CONTROL=-SIGNF( SAT, Y1+Y2*ABSF( Y2 )/(2.*SAT) ) 

C  SENSE  WHEN  OPTIMUM  CONTROL  CHANGES  SIGN 

IF  (CONTROL)  50,30,30 

C  CONTROL  POSITIVE  SET  JFLAG  FOR  LINEAR 

30  JFLAG=1 

100  FORMAT  ( 1  OH  SWITCH  AT , F 1 0. 5 , 4HFE ET ) 

PRINT  100, X( 1 ) 

C  CALCULATE  LINEAR  CONTROL 

32  CONTROL=-0.9*X( 1 ) -6 . 6 1 *X ( 2 ) -1 0 . 89*X(  3 ) -5. 81  * X ( 4 ) 

C  CHECK  FOR  SATURATION 

IF  (CONTROL-SAT)  40,40,35 
35  CONTROL=SIGNF(SAT, CONTROL) 

40  KFLAG=KFLAG+1 

IF  (KFLAG-1)  45,45,50 
45  PRINT  200, X( 1 ) 

200  FORMAT  (15H  UNSATURATED  AT  ,  F 1 0 . 5 , 4H=  E ET  ) 

50  RETURN 
ENO. 

'  END 
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SUBROUTINE  DERI  V  <X,XDOT,T) 

DIMENSION  X  (  4  )  , XDOT ( 4 ) 

C  AVERAGES  UNCOUPLED  VARIABLES  FOR  SWITCHING  LOGIC 

C  ITITIALIZE  FOR  NEW  CURVE 

IF  (T)  2,2,4 
2  JFL AG=0 
KFL AG=0 

4  XOOT ( 4 ) =- 1 .  0  *  X  ( 3 )  —  1 . 0*X { 4 ) +CONTRQL 
XDOT ( 3 ) =X ( 4 ) 

XDOT ( 2 ) =X ( 3 ) 

XDOT ( 1 )=X(2) 

SAT  =  I  00. 

C  IF  SWITCHING  POINT  PREVIOUSLY  REACHED  USE  LINEAR  CONTROL 

IF  (JFLAG)  6,6,32 

C  SWITCHING  POINT  NOT  PREVIOUSLY  REACHED 

C  CALCULATE  UNCOUPLED  STATE  VARIABLES 

6  Y 1  =X  M  )  +X  (  2  )  +  X  (  3  ) 

Y2=X(2)+X(3)+X(4) 

Y3=X ( 3 ) 

Y4=X ( 4 ) 

C  CALCULATE  OPTIMUM  CONTROL 

CONTROL=-SIGNF(SAT,Yl+Y3+(Y2*ABSF(Y2 ) +Y4#ABSF ( Y4 ) )/(2.*SAT) ) 

C  SENSE  WHEN  OPTIMUM  CONTROL  CHANGES  SIGN 

IF' (CONTROL)  50,30,30 

C  WHEN  CONTROL  CHANGES  SET  JFLAG  FOR  LINEAR 

30  JFL AG= 1 

100  FORMAT  (10H  SWITCH  AT , F 1 0 . 5 , 4HFE ET ) 

PRINT  100, X( I ) 

C  CALCULATE  LINEAR  CONTROL 

32  CONTROL=-0.9*X( 1 ) -6 . 6 1 *X ( 2 ) - 1 0 . 89*X(  3)-5.81*X(4) 

C  CHECK  FOR  SATURATION 

IF  (CONTROL-SAT)  40,40,35 
35  CONTROL=SIGNF(SAT, CONTROL) 

40  KFLAG=KFL AG+ 1 

IF  (KFLAG-1)  45,45,50 
45  PRINT  200 ,X (  1  )  ' 

200  FORMAT  (15H  UNSATURATED  AT , F 1 0 . 5 , 4H= E ET ) 

50  RETURN 
ENO 
END 
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